addpath ../Generator
addpath ../Generator/GrTheory
addpath ../PDCO
addpath ../CRNTSolvers

%Generate a random stoichiometrix matrix 
%to converge to the fixed point 
m = 100;
n = 200;
r = 3;
%Set the seed to be able to reproduce the experiment
RandStream.setDefaultStream(RandStream('mt19937ar','seed',0));
%Generate a random set of n products sets on n species, each with at most r species
Y = YGenerator(m,n,r);
%Generate the stoichiometric matrix by selecting a random set of n reactant sets where at most r species participate
S = Y-YGenerator(m,n,r);

Y = [max(-S,0),max(S,0)];
%Generate the rate constants
k = abs(randn(n*2,1));

%Generate the Ak matrix
Ak = [-eye(n),eye(n);eye(n),-eye(n)]*diag(k);


[iter,v_vecs,lmda_vecs,mass_infeas,mass_action_infeas,mu_vecs]=...
	FpIterationPDCO(Y,Ak,1e-6,5000);

fprintf('Final Mass Balance infeasibility %d\n Final Mass Action Infeasibility %d\n',mass_infeas(iter),mass_action_infeas(iter))

% plot the convergence of mu to [R F]*v
h = semilogy(1:iter, mass_infeas(1:iter),'b.-');
title('Mass balance infeasibility vs iteration count for two linkage class network')
xlabel('Iteration')
ylabel('log ||YA_k\eta||_\infty')
saveas(h,'InfeasibilityVsIterationMultiple.jpg');


